logo

This script fits state-space models to TRPD grassland bird survey data to assess species-level population trends.

Setup

Load packages & data

knitr::opts_chunk$set(echo = TRUE)
library(here)
library(png)
library(jagsUI)
library(ggplot2)
library(ggthemes)
library(tidyverse)
library(kableExtra)

source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
source(here("Scripts", "5_Prep_Species.R"))

Set parameters

# Set random seed
set.seed(42563)

# Identify species and transsects for analysis.
species_to_run <- prairie_sp
transects_to_sum <- c("CHOA01", "CHOA02", "CHOA09", "CHOA14")

# ID tag for this run and create associated output directory
tag <- "20231128_plus1_to_counts_fixplots"
mainDir <- here("Results", "991_bird_trends_species")
dir.create(file.path(mainDir, tag), showWarnings = TRUE)
outputDir <- here("Results", "991_bird_trends_species", tag)
subdirectories <- c("annual_n_estimates", "annual_n_plots", "diagnostic_plots", "species_level_estimates", "model_outputs", "synthesized_results", "Rhats", "model_parameters")
lapply(file.path(mainDir, tag, subdirectories), function(x) if (!dir.exists(x)) dir.create(x))

adjust_counts_by_1 <- TRUE

write.table(transects_to_sum, file = here("Results", "991_bird_trends_species", tag, "model_parameters", "transects_summed.txt"))

Refine species to run

transect_birds <- bird.tblBirds_Transects

# Merge BB and OA surveys that are actually the same transect
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"


# Filter bird data to only specified transects
transect_birds <- transect_birds %>%
  left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
  filter(SurveyCode %in% transects_to_sum) %>%
  mutate(
    survey_string = paste0(SurveyCode, "-", SurveyDate), # unique code for survey based on area, site #, and date
    year = sub("-.*", "", SurveyDate), # convert survey date to just year
    year = as.numeric(year)
  ) # year to numeric

# Identify days and years in which all transects in the list were actually surveyed
survey_dates <- transect_birds %>%
  group_by(SurveyDate, SurveyCode) %>%
  summarize(survey_count = length(unique(SurveyDate))) %>%
  pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  select(SurveyDate)

survey_complete_years <- transect_birds %>%
  group_by(year, SurveyCode) %>%
  summarize(survey_count = length(unique(year))) %>%
  pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  select(year)

# Re-filter bird data to only dates with complete surveys across all transects. This version forces all transects to have been run on the same day.
transect_birds <- transect_birds %>%
  filter(SurveyDate %in% survey_dates$SurveyDate)

# In years surveyed more than once, pick just the last survey from that year.
surveys_all <- transect_birds %>%
  group_by(year, SurveyCode, survey_string) %>%
  summarize() %>%
  mutate(surveycode_year = paste0(SurveyCode, "-", year))
surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
  ungroup() %>%
  select(survey_string)

# Re-filter bird data based on previous step to only include one survey day each year.
transect_birds <- transect_birds %>%
  filter(survey_string %in% surveys_selected$survey_string)

# Only include counts > 0 (remove rows where bird only detected outside survey area).
transect_birds <- transect_birds %>%
  filter(BirdCountIn > 0)

# Get species list from data and identify species observed at least 5 times.
sp_detected <- unique(transect_birds$BirdCode)
sp_detections <- transect_birds %>%
  group_by(BirdCode) %>%
  filter(BirdCode %in% prairie_sp) %>%
  summarise(
    count_indiv = sum(BirdCountIn),
    n_surveys = n(),
    n_yrs = length(unique(year))
  )
sp_detected_above_threshold <- sp_detections %>%
  filter(count_indiv >= 5) %>%
  pull(BirdCode)

# Re-filter species to run so that it only includes species actually detected in selected surveys
species_to_run <- intersect(species_to_run, sp_detected)

# Also create version with only species detected more than 5 times.
species_to_run_above_threshold <- intersect(species_to_run, sp_detected_above_threshold)

# Save these species list for use in subsequent analyses
saveRDS(sp_detected_above_threshold, here("Results", "991_bird_trends_species", tag, "sp_detected_above_threshold.RDS"))
saveRDS(species_to_run_above_threshold, here("Results", "991_bird_trends_species", tag, "species_to_run_above_threshold.RDS"))

Format data & fit model

Define function for data prep & analysis

This code creates a function for preparing time series data for species of interest in transects of interest and then fits a state-space model (assuming density independence) to the data.

# Function to prepare time series data for species and transects of interest
prep_and_fit_CHBB <- function(species_of_interest = species_to_run, transects_to_assess = transects_to_sum) {
  transect_birds <- bird.tblBirds_Transects

  # Merge BB and OA surveys that are actually the same transect
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"

  # Filter bird data to only those transects
  transect_birds <- transect_birds %>%
    left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
    filter(SurveyCode %in% transects_to_assess) %>%
    mutate(
      survey_string = paste0(SurveyCode, "-", SurveyDate), # Unique code for survey based on area, site #, and date
      year = sub("-.*", "", SurveyDate), # Convert survey date to just year
      year = as.numeric(year)
    )

  # Identify days and years in which all transects in the list were actually surveyed
  survey_dates <- transect_birds %>%
    group_by(SurveyDate, SurveyCode) %>%
    summarize(survey_count = length(unique(SurveyDate))) %>%
    pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
    ungroup() %>%
    filter(complete.cases(.)) %>%
    select(SurveyDate)

  survey_complete_years <- transect_birds %>%
    group_by(year, SurveyCode) %>%
    summarize(survey_count = length(unique(year))) %>%
    pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
    ungroup() %>%
    filter(complete.cases(.)) %>%
    select(year)

  # Re-filter bird data to only dates with complete surveys across all transects.This version forces all transects to have been run on the same day.
  transect_birds <- transect_birds %>%
    filter(SurveyDate %in% survey_dates$SurveyDate)

  # In years surveyed more than once, pick just the last survey from that year
  surveys_all <- transect_birds %>%
    group_by(year, SurveyCode, survey_string) %>%
    summarize() %>%
    mutate(surveycode_year = paste0(SurveyCode, "-", year))
  surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
    ungroup() %>%
    select(survey_string)

  # Re-filter bird data based on the previous step to only include one survey day each year
  transect_birds <- transect_birds %>%
    filter(survey_string %in% surveys_selected$survey_string)

  # Create df of years with complete data
  survey_years <- unique(transect_birds$year) %>%
    sort() %>%
    as.data.frame()
  survey_years <- rename(survey_years, year = .)

  # Count the number of species of interest across all transects
  sp_counts <- transect_birds %>%
    filter(BirdCode == species_of_interest) %>%
    group_by(year) %>%
    summarize(total_count = sum(BirdCountIn))

  # Zero fill for surveyed years
  sp_counts <- sp_counts %>%
    right_join(survey_years) %>%
    replace_na(list(total_count = 0)) %>%
    arrange(year)

  # This version of state-space model can't be fit when there are zeros in the data (populations that fit zero have no means of recovering/growing). Adding 1 to everything for now...
  sp_counts <- sp_counts %>%
    mutate(
      total_count_orig = total_count
    )

  if (adjust_counts_by_1 == TRUE) {
    sp_counts <- sp_counts %>%
      mutate(total_count = total_count + 1)
  }

  # Identify information about years for model
  yr_min <- min(sp_counts$year)
  yr_max <- max(sp_counts$year)
  year_range <- seq(yr_min, yr_max, 1) # Create a vector of survey years
  n_yrs <- length(year_range) # Counter for how many years of data
  year_range_df <- year_range %>%
    as.data.frame()
  year_range_df <- rename(year_range_df, year = `.`)

  # Add NAs for un-surveyed years
  sp_counts <- sp_counts %>%
    right_join(year_range_df) %>%
    arrange(year)

  #--------------------------

  # JAGS code for Bayesian hierarchical exponential population growth state space model
  {
    sink("Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags")
    cat("
  model {

  # Specify prior values
    init_lnN ~ dnorm(ln1, 1)
    lnN[1] <- init_lnN
    r_0 ~ dnorm(0,1)  # mean intrinsic growth rate
    sigma_r ~ dunif(0, 20) # annual SD in true pop growth rate
    obs_sd ~ dunif(0,20) # annual SD in observation (measurement) error
    # convert SD to precision (tau=1/SD^2)
    sigma_tau <- pow(sigma_r, -2)
    obs_tau <- pow(obs_sd, -2)

  # Likelihood
    for (t in 1:(n_yrs-1)){
      lnN[t+1] ~ dnorm(lnN[t] + r_0, sigma_tau)
    }

  # Observation process
    for (t in 1:n_yrs) {
      lnCount[t] ~ dnorm(lnN[t], obs_tau)
      N_est[t] <- exp(lnN[t]) # convert log(N) back to real scale
    } # close t loop

  # Derived parameters

    # Geometric mean rates of change (%/year between first and last year)
    ch <- N_est[n_yrs] / N_est[1]
    trend <- (ch^(1/(n_yrs-1)))-1

    # Trend from lognormal slope
    slope_trend <- ( n_yrs * sum(1:n_yrs * lnN) - sum(1:n_yrs) * sum(lnN) ) / (n_yrs * sum((1:n_yrs)^2) - sum(1:n_yrs)^2)

  }    # End jags model
", fill = TRUE)
    sink()
  }

  # Bundle data, identify parameters to save, submit to jags
  TRPD_data <- list(
    n_yrs = n_yrs, lnCount = log(sp_counts$total_count),
    ln1 = log(sp_counts$total_count[1])
  )
  parameters <- c("r_0", "sigma_r", "obs_sd", "N_est", "trend", "slope_trend")
  TRPD1 <- jagsUI(TRPD_data,
    inits = NULL, parameters, "Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags",
    n.adapt = 1000, n.chains = 3, n.thin = 10,
    n.iter = 300000, n.burnin = 50000, parallel = TRUE
  )

  #---------------------------------------------

  # Save model output. Models are too big to save everything, not including MCMC chains.
  mod_to_save <- TRPD1[which(sapply(TRPD1, class) != "mcmc.list")]
  saveRDS(mod_to_save, file = paste0(outputDir, "/model_outputs/", species_of_interest, "_model_outputs.RDS"))

  # Assess convergence by verifying all R-hat values < 1.1
  max.rhat <- function(x) {
    max(max_rhat <- c(sapply(x$Rhat, max, na.rm = TRUE)))
  }
  max.rhat.run <- max.rhat(TRPD1) # r.hat <1.1 suggests good convergence
  saveRDS(TRPD1$Rhat, file = paste0(outputDir, "/Rhats/", species_of_interest, "_rhats.RDS"))

  # Summarize r estimates
  r_est <- paste0("r: ", round(as.numeric(TRPD1$mean["r_0"]), 2), " (95% CI: ", round(as.numeric(TRPD1$q2.5["r_0"]), 2), "-", round(as.numeric(TRPD1$q97.5["r_0"]), 2), ")")

  # Summarize percent of r_0 posterior distribution above 0
  f_r0 <- round(as.numeric(TRPD1$f["r_0"]), 2)

  # Data frame with model estimates
  sp_values <- data.frame(
    sp = species_of_interest,
    r0_mean = TRPD1$mean$r_0,
    r0_sd = TRPD1$sd$r_0,
    r0_median = TRPD1$q50$r_0,
    r0_LCI_95 = TRPD1$q2.5$r_0,
    r0_UCI_95 = TRPD1$q97.5$r_0,
    r0_LCI_50 = TRPD1$q25$r_0,
    r0_UCI_50 = TRPD1$q75$r_0,
    f_r0 = round(as.numeric(TRPD1$f["r_0"]), 2),
    sigma_r_mean = TRPD1$mean$sigma_r,
    sigma_r_sd = TRPD1$sd$sigma_r,
    sigma_r_median = TRPD1$q50$sigma_r,
    sigma_r_LCI_95 = TRPD1$q2.5$sigma_r,
    sigma_r_UCI_95 = TRPD1$q97.5$sigma_r,
    sigma_r_LCI_50 = TRPD1$q25$sigma_r,
    sigma_r_UCI_50 = TRPD1$q75$sigma_r,
    f_sigma_r = round(as.numeric(TRPD1$f["sigma_r"]), 2),
    obs_sd_mean = TRPD1$mean$obs_sd,
    obs_sd_sd = TRPD1$sd$obs_sd,
    obs_sd_median = TRPD1$q50$obs_sd,
    obs_sd_LCI_95 = TRPD1$q2.5$obs_sd,
    obs_sd_UCI_95 = TRPD1$q97.5$obs_sd,
    obs_sd_LCI_50 = TRPD1$q25$obs_sd,
    obs_sd_UCI_50 = TRPD1$q75$obs_sd,
    f_obs_sd = round(as.numeric(TRPD1$f["obs_sd"]), 2),
    max.rhat.run = max.rhat(TRPD1),
    method = "SSMr"
  )

  write.csv(sp_values, file = paste0(outputDir, "/species_level_estimates/", species_of_interest, "_species_level_estimates.csv"))

  # Add estimates to simple dataframe
  sp_counts_w_ests <- sp_counts
  sp_counts_w_ests <- sp_counts_w_ests %>%
    mutate(
      N_est_mean = TRPD1$mean$N_est,
      N_est_sd = TRPD1$sd$N_est,
      N_est_median = TRPD1$q50$N_est,
      N_est_LCI_95 = TRPD1$q2.5$N_est,
      N_est_UCI_95 = TRPD1$q97.5$N_est,
      N_est_LCI_50 = TRPD1$q25$N_est,
      N_est_UCI_50 = TRPD1$q75$N_est,
      sp = species_of_interest,
      method = "SSMr"
    )
  write.csv(sp_counts_w_ests, file = paste0(outputDir, "/annual_n_estimates/", species_of_interest, "_annual_N_ests.csv"))

  # Plot posterior draws
  output_path <- here("Results", "991_bird_trends_species", tag, "diagnostic_plots", paste0(species_of_interest, "_diagnostic_plot.png"))
  png(output_path)
  par(mfrow = c(2, 1))

  plot(TRPD1$sims.list$obs_sd, TRPD1$sims.list$sigma_r,
    ylab = "Process error",
    xlab = "Observation error", cex = 0.3
  )
  plot(TRPD1$sims.list$r_0, TRPD1$sims.list$sigma_r,
    ylab = "Process error",
    xlab = "Process mean", cex = 0.3
  )

  dev.off()

  if (adjust_counts_by_1 == TRUE) {
    sp_counts_w_ests2 <- sp_counts_w_ests %>%
      mutate(
        N_est_LCI_95 = N_est_LCI_95 - 1,
        N_est_UCI_95 = N_est_UCI_95 - 1,
        N_est_mean = N_est_mean - 1,
        total_count = total_count - 1
      ) %>%
      mutate(
        N_est_LCI_95 = ifelse(N_est_LCI_95 < 0, 0, N_est_LCI_95),
        N_est_UCI_95 = ifelse(N_est_UCI_95 < 0, 0, N_est_UCI_95),
        N_est_mean = ifelse(N_est_mean < 0, 0, N_est_mean)
      )
  } else {
    sp_counts_w_ests2 <- sp_counts_w_ests
  }

  # Make time series plot
  plot <- ggplot(sp_counts_w_ests2, aes(x = year)) +
    geom_ribbon(aes(ymin = N_est_LCI_95, ymax = N_est_UCI_95), fill = "grey90", alpha = .6) +
    geom_line(aes(y = N_est_mean), color = "black", alpha = .7) +
    geom_point(aes(y = N_est_mean), color = "black", alpha = .7) +
    geom_point(aes(y = total_count), color = "maroon", shape = 1) +
    annotate(geom = "text", label = r_est, x = min(sp_counts_w_ests$year) + 5, y = max(sp_counts_w_ests$total_count), size = 3) +
    scale_x_continuous(breaks = seq(1960, 2021, by = 10)) +
    theme_clean() +
    labs(
      title = paste0(bird.tblBirdSpecies$BirdName[which(bird.tblBirdSpecies$BirdCode == species_of_interest)], " population dynamics"),
      subtitle = paste0("TRPD transects- ", paste(transects_to_assess, sep = "", collapse = ", "), "\nModeled estimate with 95% CI and raw counts \n", r_est),
      x = "Year",
      y = "Total count across transects"
    ) +
    theme(
      plot.title = element_text(size = 11, face = "bold"),
      plot.subtitle = element_text(size = 10, face = "italic"),
      legend.title = element_text(size = 10, face = "plain")
    )

  ggsave(plot = plot, filename = paste0(species_of_interest, "_plot.png"), device = "png", path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), width = 5, height = 3)
}

Apply function to species of interest

# # fit state space models for species of interest
# results <- lapply(species_to_run, prep_and_fit_CHBB)
# names(results) <- species_to_run
species_still_to_run <- list.files(
  path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), pattern = NULL, all.files = FALSE,
  full.names = FALSE, recursive = FALSE,
  ignore.case = FALSE, include.dirs = FALSE, no.. = FALSE
)
species_still_to_run <- gsub("\\_plot.png", "", species_still_to_run)
species_still_to_run <- unique(species_still_to_run)
species_still_to_run <- setdiff(species_to_run, species_still_to_run)

library(parallel)
numCores <- detectCores() # get the number of cores available

if (length(species_still_to_run) > 0) {
  lapply(species_still_to_run, prep_and_fit_CHBB)
  # mclapply(species_to_run, prep_and_fit_CHBB, mc.cores = numCores)
}

Visualize results

Load model results

species_level_estimates <- list.files(here("Results", "991_bird_trends_species", tag, "species_level_estimates"), pattern = "species_level_estimates.csv", full.names = TRUE) %>%
  lapply(read.csv) %>%
  bind_rows() %>%
  select(-X) %>%
  filter(sp %in% species_to_run_above_threshold) %>%
  mutate(
    r0_mean = round(r0_mean * 100, 1),
    r0_LCI_95 = round(r0_LCI_95 * 100, 1),
    r0_UCI_95 = round(r0_UCI_95 * 100, 1)
  )

Parsing annual variability

The models allow us to estimate the contribution of both observer error (undercounting/overcounting) and process error (real variability in growth rates) to the observed annual variability in populations.

error_parsing <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
  geom_linerange(aes(xmin = sigma_r_LCI_50, xmax = sigma_r_UCI_50)) +
  geom_linerange(aes(ymin = obs_sd_LCI_50, ymax = obs_sd_UCI_50)) +
  geom_label(size = 2.4, label.padding = unit(.075, "lines")) +
  theme_clean() +
  labs(
    x = "Process error (Growth variate variability)",
    y = "Observer error"
  )

error_parsing

error_parsing_repel <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
  ggrepel::geom_label_repel(size = 2.4, box.padding = 0.025, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000) +
  theme_clean() +
  labs(
    x = "Process error (Growth variate variability)",
    y = "Observer error"
  )

error_parsing_repel

ggsave(plot = error_parsing, filename = "error_parsing_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)
ggsave(plot = error_parsing_repel, filename = "error_parsing_repel_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)

Annual population indices

files <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = TRUE)
files_names <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = FALSE)
files_names <- gsub("\\_plot.png", "", files_names)

files <- files[which(files_names %in% species_to_run_above_threshold)]
files_names <- files_names[which(files_names %in% species_to_run_above_threshold)]


for (i in 1:length(files)) {
  cat(paste0("#### ", bird.tblBirdSpecies.ebird$common_name[which(bird.tblBirdSpecies.ebird$BirdCode == files_names[i])], "\n"))
  cat(paste0("![](", files[i], ")\n\n"))
  cat(paste0("**Annual trend (% change/year), 1978-2023:** ", species_level_estimates$r0_mean[which(species_level_estimates$sp == files_names[i])], "% (95% CI: ", species_level_estimates$r0_LCI_95[which(species_level_estimates$sp == files_names[i])], "% to ", species_level_estimates$r0_UCI_95[which(species_level_estimates$sp == files_names[i])], "%)", "\n\n"))
}

American Crow

Annual trend (% change/year), 1978-2023: 0.9% (95% CI: -1.5% to 3.3%)

American Goldfinch

Annual trend (% change/year), 1978-2023: 2% (95% CI: -2.6% to 7.4%)

Brown-headed Cowbird

Annual trend (% change/year), 1978-2023: -0.6% (95% CI: -5.1% to 4.6%)

Blue-winged Teal

Annual trend (% change/year), 1978-2023: -0.2% (95% CI: -11.6% to 11.5%)

Clay-colored Sparrow

Annual trend (% change/year), 1978-2023: 1.1% (95% CI: -7.6% to 10.8%)

Chimney Swift

Annual trend (% change/year), 1978-2023: 0.5% (95% CI: -1.9% to 3.3%)

Common Grackle

Annual trend (% change/year), 1978-2023: -2.4% (95% CI: -10.6% to 5.1%)

Common Yellowthroat

Annual trend (% change/year), 1978-2023: 3.8% (95% CI: -7% to 15.5%)

Dickcissel

Annual trend (% change/year), 1978-2023: 1.5% (95% CI: -1.3% to 5.1%)

Eastern Bluebird

Annual trend (% change/year), 1978-2023: 2.4% (95% CI: -3.6% to 9.1%)

Eastern Kingbird

Annual trend (% change/year), 1978-2023: 1% (95% CI: -4.8% to 6.7%)

Eastern Meadowlark

Annual trend (% change/year), 1978-2023: 6.5% (95% CI: -0.8% to 13.9%)

Field Sparrow

Annual trend (% change/year), 1978-2023: 2.3% (95% CI: -9.2% to 13%)

Grasshopper Sparrow

Annual trend (% change/year), 1978-2023: 0.8% (95% CI: -4.7% to 6.9%)

Henslow’s Sparrow

Annual trend (% change/year), 1978-2023: 5.2% (95% CI: -6.8% to 16.9%)

Killdeer

Annual trend (% change/year), 1978-2023: 1% (95% CI: -3.6% to 5.8%)

Mourning Dove

Annual trend (% change/year), 1978-2023: -2.3% (95% CI: -10% to 5.6%)

Northern Cardinal

Annual trend (% change/year), 1978-2023: 0.5% (95% CI: -5.4% to 6.2%)

Ring-necked Pheasant

Annual trend (% change/year), 1978-2023: 0.9% (95% CI: -1.5% to 3.9%)

Red-winged Blackbird

Annual trend (% change/year), 1978-2023: 1.2% (95% CI: -11.2% to 13.9%)

Savannah Sparrow

Annual trend (% change/year), 1978-2023: -1.2% (95% CI: -7.7% to 5.1%)

Sedge Wren

Annual trend (% change/year), 1978-2023: 2.7% (95% CI: -5.5% to 11.6%)

Song Sparrow

Annual trend (% change/year), 1978-2023: -1.4% (95% CI: -13.9% to 10.9%)

Tree Swallow

Annual trend (% change/year), 1978-2023: 5.7% (95% CI: -3.9% to 15.4%)

Vesper Sparrow

Annual trend (% change/year), 1978-2023: -2.2% (95% CI: -9.1% to 6.5%)

Western Meadowlark

Annual trend (% change/year), 1978-2023: -1.7% (95% CI: -9.6% to 5.5%)

Willow Flycatcher

Annual trend (% change/year), 1978-2023: 0.3% (95% CI: -4.6% to 4.8%)


Print session info:

sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] auk_0.7.0              jagsUI_1.5.2           png_0.1-8              readxl_1.4.3          
##  [5] kableExtra_1.3.4       ggthemes_5.0.0         transformr_0.1.3       gifski_1.12.0-2       
##  [9] gganimate_1.0.8        ggspatial_1.1.9        ggrepel_0.9.4          rgdal_1.6-4           
## [13] sp_2.1-2               RODBC_1.3-23           RVAideMemoire_0.9-83-7 pairwiseAdonis_0.4.1  
## [17] cluster_2.1.4          goeveg_0.6.5           vegan_2.6-4            lattice_0.20-45       
## [21] permute_0.9-7          sf_1.0-14              lubridate_1.9.3        forcats_1.0.0         
## [25] stringr_1.5.1          dplyr_1.1.4            purrr_1.0.2            readr_2.1.4           
## [29] tidyr_1.3.0            tibble_3.2.1           ggplot2_3.4.4          tidyverse_2.0.0       
## [33] here_1.0.1            
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.1-0   class_7.3-20       rprojroot_2.0.4    htmlTable_2.4.2    base64enc_0.1-3   
##  [6] rstudioapi_0.15.0  proxy_0.4-27       farver_2.1.1       bit64_4.0.5        fansi_1.0.5       
## [11] xml2_1.3.5         splines_4.2.2      cachem_1.0.8       knitr_1.45         spam_2.10-0       
## [16] Formula_1.2-5      jsonlite_1.8.7     rjags_4-15         compiler_4.2.2     httr_1.4.7        
## [21] backports_1.4.1    assertthat_0.2.1   Matrix_1.6-4       fastmap_1.1.1      cli_3.6.1         
## [26] tweenr_2.0.2       htmltools_0.5.7    prettyunits_1.2.0  tools_4.2.2        dotCall64_1.1-1   
## [31] coda_0.19-4        gtable_0.3.4       glue_1.6.2         maps_3.4.1.1       Rcpp_1.0.11       
## [36] cellranger_1.1.0   jquerylib_0.1.4    vctrs_0.6.5        svglite_2.1.2      nlme_3.1-160      
## [41] xfun_0.41          rvest_1.0.3        timechange_0.2.0   lpSolve_5.6.19     lifecycle_1.0.4   
## [46] MASS_7.3-58.1      scales_1.3.0       vroom_1.6.4        ragg_1.2.6         hms_1.1.3         
## [51] fields_15.2        yaml_2.3.7         gridExtra_2.3      sass_0.4.7         rpart_4.1.19      
## [56] stringi_1.8.2      highr_0.10         e1071_1.7-13       checkmate_2.3.0    epuRate_0.1       
## [61] systemfonts_1.0.5  rlang_1.1.2        pkgconfig_2.0.3    evaluate_0.23      htmlwidgets_1.6.4 
## [66] labeling_0.4.3     bit_4.0.5          tidyselect_1.2.0   magrittr_2.0.3     R6_2.5.1          
## [71] generics_0.1.3     Hmisc_5.1-1        DBI_1.1.3          pillar_1.9.0       foreign_0.8-83    
## [76] withr_2.5.2        mgcv_1.8-41        units_0.8-5        nnet_7.3-18        crayon_1.5.2      
## [81] KernSmooth_2.23-20 utf8_1.2.4         tzdb_0.4.0         rmarkdown_2.25     progress_1.2.3    
## [86] grid_4.2.2         data.table_1.14.8  webshot_0.5.5      digest_0.6.33      classInt_0.4-10   
## [91] textshaping_0.3.7  munsell_0.5.0      viridisLite_0.4.2  bslib_0.6.1
 




By Sam Safran